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ABSTRACl 


Nuincrica!  computations  have  been  performed  to 
investigate  the  interaction  of  two  counter-rotating  vortices  with 
a  free  surface  in  a  viscous,  incompressible  fluid  where  the 
motion  is  considered  two-dimensional  and  laminar.  The 
decaying  vortices  deform  the  free  surface  during  their 
approach  and  they  can  display  the  phenomenon  of 
“rebounding."  The  numerical  scheme  is  based  on  the 
Navier-Stokes  equations  and  uses  boundary-fitted  coordinates 
to  accommodate  the  locally  high  vorticitv  of  the  moving 
vortices  and  the  nonlinear  deformation  of  the  free  surface. 

Al)  M I M  SIR  ATI  V  i:  I  MO  RM  A  HON 

I'his  project  was  supported  by  the  Office  of  the  C'hief  of  Naval 
Research,  and  administered  by  Dr.  I-dwin  P.  Rood,  Fluid  Dynamics  Prostram 
( 1 132F)  under  ONR  Contract  No.  N0001489-WX-24020. 

INTRODUCTION 

I'vvo  counter-rotating  vortices  in  a  two-dimensional  viscous  fluid 
approach  through  self-induction  a  free  surface  and  interact  with  it.  The  fluid  is 
considered  incompressible  and  the  laminar  flow  transient,  with  the  fluid 
ultimately  coming  to  rest.  T'his  scenario  may  be  envisioned  as  an 
approximation  to  a  pair  of  vortex  filaments  in  the  far  wake  of  a  slowly  moving 
body  under  water. 

Various  models  of  different  relevance  to  this  problem  can  be  found  in 
the  literature.  The  classical  Lamb  solution  describes  the  potential  flow 
generated  by  a  pair  of  counter-rotating  point  vortices  approaching  or  leaving  a 
straight  boundary.  The  paths  of  the  vortices  are  given  by 

.v  -f  y“  =  4.vy‘-  (1) 

in  a  Cartesian  coordinate  system  (.v,  y)  with  the  unit  distance  between  the  two 
vortices  far  away  from  the  boundary.'  The  behavior  of  a  pair  of  vortices  with 
an  elliptic  core  was  studied  by  Saffman.“ 

In  a  viscous  fluid  an  additional  boundary  condition  must  be  specified  for 
the  tangential  velocity  component.  Under  nonslip,  vorticity  of  opposite  sign  to 
the  approaching  vortices  will  be  generated,  and  the  paths  wilt  be  different  from 
that  described  by  FT].  (1).  Flow  separation  can  occur  at  the  boundary  with  the 
subsetjuent  development  of  secondary  vortices.  The  primary  vortices  will  turn 
away  from  the  boundary,  a  phenomenon  that  is  called  “rebounding''.^  ''  A 
numerical  computation  of  the  boundary  layer  at  the  wall  was  performed  by 
F'Tsoy  and  Walker.^ 

Rebounding  was  also  observed  experimentally  by  Barker  and  Crow''  at  a 
water  surface.  No  information,  however,  was  given  on  the  shape  of  the 
surface.  The  vortices  were  generated  by  a  moving  and  then  abruptly  retracted 
plate  in  a  water  tank.  Different  methods  of  vortex  generation  were  used  by 
Sarpkaya  and  I Ientlers(m'’  and  by  Willmarth  et  al.^  '  no  stuilied  the 
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disturbance  ol  llic  water  surtacc  in  three  dimensions  by  the  steady  movement 
of  a  liydrofoil  and  by  counter-rotating  Haps,  respectively. 

/\  numerical  approach  was  used  by  Sarpkaya  et  al.,‘^  I'elste.'^ 
Marcus,''*  and  Willmarth  et  al.^  for  two-dimensional  potential  How  with  twt) 
counter-rotating  vortices  approaching  a  nonlinear  free  surface. 

Ihe  corresponding  viscous  problem  was  tackled  by  Peace  aiul  Riley" 
for  a  plane  surface  with  either  a  noiislip  or  a  perfect-slip  condition.  I'or  the 
initial  Ilow  development  they  used  a  viscous  inner  solution  and  an  inviscid 
outer  solutivut.  After  this  initial  phase,  a  fmitc-diffcrcncc  scheme  for  the 
Na\ier-Su'ke>  equations  was  tipplied. 

In  this  report  the  viseous-llow  problem  is  solved  for  »  nonlinear  free 
surface  whieli  produces  vorticity  according  to  the  formula:  Vorticity  ec|uals  twc) 
times  curvature  times  tangential  velocity  at  the  free  surface.  y\  linite- 
tiillerenee  scheme  is  used  lor  the  entire  time  span  during  which  the  vortices 
approach  the  free  surface. 

s  i  ah  .Mi;.M  ov  i  he  problem 

.\  pair  of  vortices  of  equal  strength  t:  but  opposite  sign,  a  distance  a 
apart,  approaches  through  self-induction  an  initially  undisturbed  free  surfaee  at 
y'  0  in  the  C  artesian  coordinate  system  (.v',  y')  (Fig.  1).  With  the 
correspomling  velocity  components  u'  and  v'  the  time-dependent  How  field  of 
the  incompressible  Newtonian  tluid  is  described  by 


■  ni' 

,ihi 

+  u  — r  -r  '■ 

Du' 

1  Dp' 

^  D-u' 
+  'A 

+  ^)  > 

(2) 

.'y' 

d.v 

Dy ' 

l>  Dx' 

Dx'' 

dy’' 

.  y’ 

.  Dy' 

1  Dp' 

+  a 

+  T  )  ' 

(3) 

lit' 

<)x' 

Dy ' 

p  <^y' 

Dx'- 

dy'- 

Du' 

Dv' 

=  0  , 

(4) 

Dx ' 

Dv ' 

with  p\  p.  and  g  the  time,  total  pressure,  density,  kinematic  viscosity, 
anti  eamstant  of  gravity,  respectively. 

The  free  surface  is  described  by  y' =  Y'(x',  (')  and  is  part  of  the 
solution.  On  this  free  surface  the  boundary  eonditions  are 

... 


(6) 


,  ,  „  du\dY'  ,  ,  Ou'  ,  Ov\  „ 

^  +  "‘Tv  +  7v>  =  "  ■ 


/  ,  ^  d\'\  ,  /dll'  dv'.dY'  . 

-  -"Tv'  +  "‘Tv  +  Tv >Tv  =  "  ■ 


Surface  tension  is  neglected.  The  other  boundary  conditions  are 
y' =  —CO  :  u' =  c' =  0  ,  p' =  +00  , 

.v'  =  0  ,  -,v  <  v'  <  r'  :  =  0  ,  n'  0  .  ^  =  0  . 

'  d.V  d.V 

.v'  =  v-cc  ,  — ec  <  v'  <  0  :  u'  —  e’  =  0  ,  />'  =  —  /igy'  . 


(7) 

(«) 

(9) 

(10) 
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Since  tlie  How  field  is  symmetric  about  the  line  .v’  —  0,  only  the  half 
plane  .v’  >  0  with  one  vortex  at  .v'  =  .vV  and  y'  =  y'^,  is  considered  (I'ig.  1). 

I'he  following  assumptions  are  made  for  the  initial  conditions  at  t'  —  I'o: 

1.  I  he  free  surface  is  undisturbed. 
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A  point  vortex  with  strength  k  at  the  position  .v'  =  x’^  and  y'  =  vV  is 
introduced  whose  flow  field  is  irrotational  and  is  described  in  llie 
quadrant  .v'  >  0,  y'  <  0  by 


d'  +  /(/■'  = 


-ih-  log 


(t:' —  +z'v) 

(z'-Fixz’+r:) 


(11) 


where  •'  =  .v'  +iy'  and  <;)'  and  V’’  are  the  potential  function  and  the  stream 
function,  respectively. 

3.  I'he  analytic  solution  (11)  of  the  How  field  does  not  satisfy  the  boundary 
condition  (7).  This  inconsistency  has  only  a  small  numerical  effect  on 
the  first  few  time  steps. 


4.  Since  the  point  vortex  represents  a  singularity  with  infinite  vorticity  in  the 
flow  field,  it  is  practical  to  assume  a  certain  decay  of  the  point  vortex  to 
accommodate  it  for  a  finite-difference  grid.  If  this  decay  is  confined  to  a 
local  area  about  the  center  of  the  vortex  so  that  there  is  no  appreciable 
effect  on  the  second  vortex,  that  is,  on  the  symmetry  line  (or  on  the  free 
surface),  the  Hamcl-Oseen  solution  can  be  applied:’*  '" 

0 


-7l' 


with  I 


f 


exp{- 

^  r 

=  ^ ir  -y  r"  and  r'^  =  (.v'  - 


- - ^ - )1 

•v'v)'  +(v’-yV)'- 


(12) 


It  is  convenient  to  introduce  the  following  dimensionless  quantities: 


3 


(1.') 


:.v’.  v',  V)  =-  aix.  V.  Y)  ,  t'  =  ,  (//',  y')  Vn{u,  r) 

>'  0 

---  I>\':){P  -y  /I- n)  .  ^  ~a 


with  //  the  initial  depth  ot'  the  vortices,  V'd  the  initial  translational  velocity  ol 
the  vortex  pair  (it  may  be  mentioned  that  this  scaling  makes  the  unit  distance 
in  b.q.  (1)  only  approximate,  that  is,  0.99),  and  P  the  scaled  dynamic  pressure. 
The  ilirnensionless  How  parameters  are  the  Fr(.>iide  and  the  Reynolds  numbers 
I  '-  and  Iii\  respectively; 


Y'i.a  K 

r  i' 


(  14) 


.Moreover,  the  inlinite  domain  of  integration  is  approximated  by  a  finite 
ilomain  lor  numerical  reasons.  Also,  the  inertial  terms  are  expresserl  in 
eonseiwation  lorm.  Then,  the  initial-boundary  value  problem  is  delineil  by 
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(17) 

dx 

dy 

with  the  following  boundary  conditions; 


y  -  y; 
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y 
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dx 
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(19) 

(20) 


xj  ;  p  =  []  .  u,  r ;  obtained  by  2nd^order  extrapolatimi 
along  coordinate  line  into  the  interior  , 

df  „  <)P 


X  --  0  ,  — v/  <  y  Y  y  ■ 


0  ,  U  :  0 


0  , 


Hx  <>x 

-^Xf  ,  — V/  <  \'  '  0  ;  P  ^  0  ,  n,  r;  obtained  by  2nd— ortler 
extrapolation  along  cooialinate  line  inter  thi.  interior  . 


(21) 

(22) 

(2.V) 


1  he  initial  comlitions  at  I  t;  are; 
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Position  of  vortex:  a\.  ==  1/2,  \v  =  — 

I'he  flow  (icid  is  irrolational  and,  according  to  l  iq.  (11), 


d+i(l’=  —/log 


(-~  ~  ~v)("  +  ^v) 

(c  ^  +  Cv) 


(24) 


except  for  the  vicinity  of  the  vortex  center  witii  the  circular  boundary  r/  : 


expi¬ 


re  Re 


4U/,  -  /,)) 


n. 


(25) 


I'or  a  certain  r/  ,  which  is  determined  by  the  permissible  error,  f/  can 
be  computed  and  is  for  the  examples  in  this  report  //  —  0.125  with  /q  = 
This  //  is  then  chosen  as  the  start  of  the  computation:  t  ~  ti  —  0  with 
to  =  -0.125. 


NUMERICAL  TECHNIQUE 

I'he  numerical  solution  of  the  initial-boundary  value  problem,  as 
expressed  in  Pqs.  (15)  through  (25),  is  carried  out  with  the  aid  of  a  finite- 
difference  technique  and  boundary-fitted  coordinates. 

Boundary-fitted  coordinate  transformation 

Fig.  lb  is  a  schematic  drawing  of  how  the  physical  s'pace  is  mapped  onto 
the  computational  domain.  Only  the  coordinate  lines  which  form  the 
boundaries  of  the  two  regions  are  drawn.  The  coordinate  lines  in  physical 
space  are  mapped  onto  a  uniformly  spaced  Cartesian  mesh  with  a  unit  mesh 
spacing  in  each  coordinate  direction. 

As  the  flow  field  evolves  in  time,  the  grid  in  physical  space  will  move 
and  its  coordinate  lines  will  be  attracted  to  regions  of  high  flow  gradients 
through  the  use  of  an  adaptive-grid  technique  to  be  discussed  later.  However, 
the  ('artesian  grid  in  computational  space  always  remains  fixed  and  unif  rm. 
This  is  the  major  advantage  of  using  a  mapping. 

For  Re  =  10,  f’r  =  1.125,  the  physical  region  extends  from  .v  =  0  to 
n.O  and  from  y  =  —6.0  to  the  free  surface  which  is  initially  at  y  =  0.  This 
physical  region  is  mapped  onto  a  computational  space  with  a  Cartesian  grid 
consisting  of  91  equally-spaced  points  in  both  the  and  //-  coordinate 
directions.  For  the  two  Re  =  50  cases,  the  physical  regions  extend  from 
.V  =  0  to  10.8  and  from  y  =  —6.0  to  the  free  surface.  These  physical  regions 
are  mapped  onto  a  compulatioftal  space  with  a  Cartesian  grid  consisting  of  157 
points  in  the  Cvjordinate  direction  and  135  points  in  the  //—coordinate 
direction.  Figs.  2  through  4  display  representative  physical-space  grids  at 
various  times. 


I'hc  curvilinear  coordinates  (^,  r/)  are  obtained  as  solutions  of  the  two 
elliptic  partial  differential  equations  with  the  physical-space  coordinates  (.v,  y) 
as  independent  variables 


+sSy  =  (d  +  s^?)  n)  , 

'/.«  +  f/vy  =  i'll  +  'll)  Q*(^,  '/)  • 


(26) 


(27) 


However,  since  all  calculations  are  to  be  done  in  the  rectangular  computational 
domain,  these  two  elliptic  partial  tliffcrential  equations  are  transformed  by 
interchanging  the  dependent  and  independent  variables.  As  a  result,  the 
physical-space  coordinates  (.v,  y)  are  solved  in  terms  of  the  computational 
space  coordinates  (^,  //)  at  each  time  step.'’  The  transformation  yields,  with 
partial  derivatives  now  symbolized  by  subscripts  for  simplicity. 


O.V^^  —  2d.V£.,  -h  +  o/’-’-.Vc  -p  ''iQ'X,.  —  0  , 

oy.-;  —  2.(v,-,,  -f  ~,y,„,  -f-  nP'-y.-  -P  yQ'-'y,!  =  0  , 


(28) 

(29) 


with 


A  =  .v-  -py- 


.1=  X:X,I  yyyy,,  (30) 

■>  *> 

"/  =  +  yi  ■ 

1  \)r  later  use,  the  Jacobian  of  the  transformation  is  given  by 

y  =  .v^y„  -  . 

riie  form  of  the  control  functions  P*  and  Q*  for  the  coordinate  system 
will  be  presented  when  adaptive-grid  generation  is  discussed. 

The  Navicr-vStokes  equations  (15)  and  (16)  are  in  curvilinear  coordinates 
(s'. 

u,  -  .v,(y„a£  -y^u„)/J  -  y,(.v^w,,  -x,,u{)/J 


+  l.'^,(t/^)£  -  .Vf(w')„l/y  -P  [.v^(//i’)„  -x,,{uv)(]/J  -f  (y„P^  ~y(P„)/J 
---  i<Hnf  +  yu,,,,  +  a*u,i  +  T*ii^)/ReJ'  ,  (31) 

'V  ~x,(y„vc  -y^i'„)/y  -yiix^y,,  -x„\’^)/J 
-Ply„(m’)^  -yiiuv)„]/J  +[x^(v~)„  -x„{\'^)^\/J  -f  (.v^P„  -x„P^)/J 
=  (nv^^  -  -P  yv,,,,  -f  rrH’ ,,  -p  rH-^)/ReJ~  ,  (32) 


where 


6 


(33) 


//)  .  r*  -  //) 


The  timc-ckMivativcs  have  also  been  transbirmccl  in  tliese  equations. 
I'hus,  time  derivatives  in  I'qs.  (.^1)  and  (32)  arc  taken  with  ^  and  ;/  fixed,  while 
those  in  I  'qs.  (15)  and  (!'^}  v  ere  taken  will*  .i  and  y  fixed.  I  his  tninsfonnution 
of  time-rierivatives  allows  the  computation  to  be  dime  on  a  fixeti  grid  in  the 
transformed  plane  even  thi>ugh  the  physical  grid  is  in  motion  ilue  to  the 
movement  ol  llie  free  surface.  I'he  terms  involving  .i,  and  \',  in  liiis.  (31)  and 
(32)  hence  occur  because  of  the  miwing  physical-space  grid.  This  procedure 
was  adopted  from  Shanks.'* 


The  continuity  et|uation  (17)  is  replaceil  by  an  equation  with  pseudo¬ 
compressibility  for  numerically  conserving  mass  at  each  physical  time  step; 


,>r 


.'It 


+  3^- 


iky 

tx 


0  . 


(34) 


In  this  method  pseiulo-time  steps  for  pscudi*-limc  r  are  rcquireil  to 
satisfy  bq.  (17)  at  each  physical  time  step  This  approach  can  be  vieweti  as 
an  iteration  procedLire  in  pseudo-time  7  to  calculate  each  physical  time  step  A/. 
The  method  has  been  successfully  applied  recently  by  several  researchers*'  ’'’ 
to  conqnite  the  two-dimensional,  incompressible  Navier-.k'.okes  equations  for 
time-dependent  flows.  'I'liis  technique  was  first  introduced  by  ('horin*^  for 
obtaining  solutions  to  the  steady-state,  incompressible  Navier-Stokes  equations 
and  is  characteri/ed  by  the  existence  of  moving  pseudo-pressure  waves.  They 
die  out  in  pseudo-time  leaving  a  divergence-free  velocity  field  at  steady  state. 
The  method  of  employing  artificial  compressibility  has  been  used  successfully 
by  many  authors  for  computing  steady-state  incompressible  Navier-Stokes 
solutions  (see.  for  example,  Kwak  et  al.”^). 


'I  he  parameter  J*  can  assume  values  between  zero  and  ten  but  is  usually 
chosen  to  be  one,  and  so  it  is  in  the  present  work.  Ihen,  liq.  (34)  is  in  the 
transformed  computational  space 

OF  1 

—  +  7Cv,;W?  -.v?w„  +x^y,,  -  .v„f^)  =  0  .  (35) 

(JT  J 


fire  movement  of  the  grid  during  the  p.scudo-time  step  Ar  is  so  small 
that  the  terms  involving  .v^  and  y.  have  been  neglected  in  f,q.  (35). 

The  kinematic  boundary  condition  (18)  at  the  free  surface  is  for 
coordinates  in  computational  space 

OY  - 

V  =  f :  m 

Ot  .V 


Ihe  subscript  ,v  refers  to  the  physical  space  grid  points  at  the  free 
surface  that  are  nirt  pernritted  to  move  in  the  .v  -direction.  An  equivalent 
kincin.iiic  cojulilion  to  liq.  (36)  is  given  in  the  following  equation  in  which  the 
grid  )roints  at  the  free  surface  are  allirwed  to  move  in  the  .v —direction; 


r' 


y; 


HY 

Ot 


i’ 


d.V 

- =  n 

in 


(37) 
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Both  equivalent  kinematic  conditions  (36)  and  (37),  of  course,  permit 
the  "rid  points  at  the  free  surface  in  physical  space  to  move  in  the 
\'  —direction. 

I'he  free  surface  in  physical  space  maps  onto  a  constant  //—line  in 
compulational  space,  as  seen  in  I’ig.  lb.  h'or  a  constant  //—line,  the  free 
surface  conditions  (19)  and  (20),  written  in  computational-space  coordinates, 
are 


y  =  T: 

=  -1.-)//.  -Av  -ReJy.iP  - 
/ 

7^)1  , 
ir~ 

(38) 

=  ~1A/,-  -I-  +  ReJxc{P  — 

1 

• 

(39) 

In  addition  to  the  kinematic  condition  for  Y  at  the  free  surface,  three 
conditions  are  needed  to  determine  P,  u,  and  r  at  the  free  surface.  Condition 
(38)  is  retained  to  compute  u  whereas  a  linear  combination  of  liqs.  (.38)  and 
(.39)  is  derived  that  yields  the  condition  for  /’.  The  velocity  component  i’  is 
then  obtained  from  the  continuity  equation  to  conserve  mass  at  the  free 
surface.  I  he  conditions  for  C  and  f  at  the  free  surface  v  =  Y  are 

y„U(  -y'iU,,  +->-T'9/  -  -^'/'T  =  ^  •  (4^) 

At  .V  =  0,  the  symmetry  conditions  (22),  after  being  written  in 
computational-space  coordinates,  are  used  to  compute  v  from  the  Navier- 
Stokes  equation  (32)  and  P  from  Cq.  (35).  'I'he  y -coordinate  at  .v  =  0  is 
obtained  from  liq.  (29).  Of  course,  u  and  x  are  always  zero  on  the  symmetry 
line  .V'  =  0,  including  the  symmetry  point  at  the  free  surface.  At  this  special 
point,  Y,  P,  and  f  are  obtained  from  Eqs.  (36),  (40),  and  (41),  respectively, 
with  the  symmetry  conditions  (22)  built  in. 

Fhe  boundary  conditions  (21)  and  (23)  are  used  with  second-order 
extrapolations  along  interior  coordinate  lines  in  computational  space  to  obtain 
u  and  V  at  these  boundaries.  For  instance,  at  y  =  —yi,  u  is  obtained  from 
values  of  u  at  the  two  grid  points  on  a  constant  line  closest  to  the  boundary. 
Similar  considerations  apply  to  v  and  to  the  boundary  of  Eq.  (23),  except  that 
constant  //—lines  are  used  there. 

Adaptive  gridding  is  used  in  this  report  by  giving  a  special  form  to  the 
coordinate  system  control  functions  P*.  Q*,  which  appear  in  Eqs.  (28),  (29), 
(31),  and  (32).  The  basic  idea  is  to  use  the  equi-distribution  of  a  weight 
function  along  arclength  elements  in  the  physical-space  grid.'^  'I’hese  equi- 
distribution  laws  for  weight  functions  ic,  and  u't  along  arcle/igth  elements  on 
co/istant  //—  and  lines,  respectively,  are 


1 


(.17  -t-  yj )  “  »’  I  =  const  , 

(42) 

X 

(.\7,  -f  y:,) '  ux  =  const  . 

(4.3) 
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Weight  functions  are  usually  taken  to  be  functions  of  the  How  gradient, 
and  they  are  chosen  here  to  be 


U'l  =  (1  +/i- 


(1  +qW- 


Wl+AVi 


(44) 


ii'i 


(1  +/?-| 


'  iin 


i)Vr 


A-qj, 


(45) 


where  q 


v7 


+  r- 


One  observes  from  Eqs.  (42)  through  (45)  that  the  spacing  of  the 
arclength  between  grid  points  in  physical  space  will  be  small  if  the  gradient  of 
the  local  flow  speed  q  is  high.  The  grid  adapts  to  the  locally  high  flow  gradient. 
If  the  local  flow  gradient  is  zero,  uniform  spacing  will  occur. 

Anderson''^  has  shown  that  if  P*  and  Q*  have  the  form 


Q*  = 


(»*’2)„ 

U’2 


(46) 


the  mesh-generating  equations  (28)  and  (29)  approximate  the  equi-distribution 
laws  (42)  and  (43),  respectively.  In  the  computations,  coefficients  Cj  and  CA 
are  actually  added  to  the  right  sides  of  Eq.  (46).  C]  and  Cz,  which  are 
constant  in  time  but  vary  spatially,  are  the  initial  P*  and  Q*,  respectively,  for 
the  Initial,  non-uniform  Cartesian  grid  in  physical  space  obtained  from  the 
INMESH  program. This  grid  is  used  at  the  very  start  of  the  flow 
computation  at  /  =  0  with  adaptive  gridding  then  applied  immediately. 

On  the  symmetry  line  x  =  0,  the  weight  function  Wz  is  given  in  the 
simpler  form 

u.'2  =  1  +  A\’l  .  (47) 

Only  Q*  in  Eq.  (46)  and  Eq.  (47)  arc  required  in  conjunction  with  Eqs. 
(29)  and  (32)  to  compute  y  and  v,  respectively,  on  the  symmetry  line  .v  =  (). 

Except  for  the  symmetry  line,  adaptive  gridding  is  used  only  in  the 
interior  of  the  flow  region. 


Finite-Difference  Technique 

All  spatial  derivatives,  including  one-sided  derivatives  at  the  boundaries, 
are  replaced  by  finite-difference  operators  of  second  order  in  the 
computational  space.  Thus,  it  remains  to  discuss  only  the  implicit  time- 
differencing  procedure  for  the  initial-boundary  value  problem  that  consists  of 
Eqs.  (21)  through  (25),  (28)  through  (33),  (35),  (36)  or  (37),  (3«),  (40),  (41). 
and  (46).  I’he  weight  functions  in  E!q.  (46)  are  provided  by  Eqs.  (44),  (45)  or 
(47). 

The  dynamic  pressure  field  E  at  /  =  t/  =  0  is  obtained  by  solving  a 
Poisson  equation  for  P  in  terms  of  the  initial  velocity  field,  Eqs.  (24)  and  (25). 
I’his  is  the  only  time  a  Poisson  equation  for  P  is  used. 
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The  lollowiiig  equations  represent  the  implicit  time  diirerencing 
procedure  lor  advancing  the  liow  solution  {l\  u,  c,  .v,  y)  in  the  interior  lor  the 


physical  time  step  _V  =  from 

using  pseudo-time  steps  Ar  =  —  r"' 

physical  time  level  n  lo  level 

n  -i  1  by 

y/l+l.  f/l+l  _  ^  ^  +  1  ,  »)  +  1 

^./i  +  l.  m+1  p^.n  Q*"^ 

(48) 

+  ^.)l  + 1 .  m  +  1 

y''  +  ''"'  +  i,  />*'■,  CP-")  . 

(49) 

pii  +  l.  m  +  1  _  +  in 

Ar 

V'T*)"  ^ ' 

(50) 

+  m  +  1  _ 

t- 1 .  in  4  1 

(51) 

A/ 

-  1 

|,;i  +  1 .  (H  + 1  _  |,'i 

T  1 .  m  +  1 

(52) 

A/ 

=  ^2 

liqs.  (48)  through  (52)  are  obtained  I'rom  Ixis.  (28),  (29),  (35),  (31),  and 
(32),  respectively.  Superscripts  refer  to  time  levels.  I'he  nonlinear  functions 
y>i  and  ^2  I'Ms.  (48)  and  (49)  are  functions  of  the  latest  updated  values 
^.«  +  i.m  +  i^  ^,«  +  i.m+i  .^  neighboring  points  after  using  spatial  differencing. 
In  I'.qs.  (51)  and  (52)  the  time  derivatives  for  u  and  i'  are  on  the  left  side  of  the 
equations  while  all  remaining  terms  of  liqs.  (31)  and  (32)  (in  difference  form) 
are  ittcluded  in  the  functions  r\  and  on  the  right  side  of  Hqs.  (51)  and  (52). 
Physical  time  derivatives  occur  in  liqs.  (31),  (32),  (36)  or  (37)  and  are 
differenced  to  first  order  according  to 


I 


/I  + 1 .  m  + 1 


-A/ 


(53) 


where  /  stands  for  u,  v,  .v,  y,  or  V.  All  spatial  derivatives  of  the  initial¬ 
boundary  value  problem,  after  second-order  finite  differencing,  are  evaluated 
by  using  the  latest  available  updated  values  of  the  implicit  scheme. 


'I'he  flow  solution  at  the  new  physical  time  level  /i  -t-  1  is  obtained  wh.en 
the  convergence  criteria 


1  + 1 ,  m  + 1  1 

>  <1  - 

(54) 

r/i  + 1 ,  m  -f  1  r»i  +  U  m 

iJ - L - 1 

'  y'M-fl.  m  ' 

<  <2 

(55) 

arc  satisfied  for  f  =  u,  v,  .v,  and  v  at 

all  grid  points  of 

the  computational 

space.  'Then,  +  l„  p^js 

.  (54)  and  (55)  rj 

1  and  <2  are  small 

specified  parameters. 

Rogers  and  Kwak'^  used  the  implicit  scheme  of  Tqs.  (50)  through  (52) 
for  How  problems  with  fixed  geometry.  Rqs.  (51)  and  (52)  can  be  viewed  as  an 
Tiller  backward  scheme  for  the  Navicr-Stokes  equations.  The  convergence  of 
P  in  pseudo-time  ensures  conservation  of  mass  at  each  physical  time  step  as 
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discussed  earlier.  This  process  is  also  seen  from  Fq.  (50). 

Rogers  and  Kwak^^  used  upwind  differencing  for  the  convective  terms 
of  F.qs.  (51)  and  (52)  which  is  necessary  for  high  Reynolds  numbers.  In  this 
work  a  central  difference  operator  was  applied. 

A  “four-color”  scheme  (Fig.  Ic)  is  used  in  the  interior  of  the 
computational  space  for  Eqs.  (48)  and  (49)  and  for  Fqs.  (50)  through  (52)  in 
obtaining  the  latest  values  with  superscripts  n  +  I,  m  +  \  (updates)  for  x,  v 
and  P,  u,  and  c,  respectively.  'I’he  use  of  such  a  scheme,  which  can  be 
vectorized,  resulted  in  an  order  of  magnitude  increase  in  computer  speed  on 
the  ('ray-XMP  24  on  which  the  computations  were  performed. 

The  “four-color”  scheme,  as  applied  to  liqs.  (50)  through  (52),  for 
example,  consists  of  obtaining  updates  for  P,  u,  v  simultaneously  at  all  the  o 
points,  then  at  all  the  points,  the  x  points,  and  the  A  points,  in  that  order. 
The  latest  available  updates  are  used  in  this  process. 

The  computational  cycle  for  one  complete  pseudo-time  step  iteration 
consists  of  (a)  applying  the  “four-color”  scheme  to  Fqs.  (48)  and  (49)  followed 
by  obtaining  the  latest  updates  for  v  at  successive  points  along  the  symmetry 
boundary  Fq.  (22)  from  Fq.  (49);  (b)  applying  the  “four-color”  scheme  to 
Eqs.  (50)  through  (52);  (c)  obtaining  updates  for  P  and  e  at  successive  points 
along  the  symmetry  boundary  Fq.  (22)  from  Fqs.  (50)  and  (52),  respectively; 
(d)  obtaining  updates  for  P,  u,  v,  and  Y  at  successive  points  along  the  free 
surface  from  Fqs.  (40),  (38),  (41),  and  (36),  respectively;  and  (e)  obtaining 
updates  for  u  and  f  at  successive  points,  first  along  the  boundary  .v  =  .V/,  from 
Fq.  (23)  and  then  along  the  boundary  _v  =  — V/,  from  Fq.  (21). 

At  the  completion  of  this  computational  cycle,  after  the  latest  updates 
for  .V',  y,  u,  and  c  satisfy  the  convergence  criteria  of  Eqs.  (54)  and  (55)  at  all 
points,  these  updates  are  the  solution  at  the  new  time  level  n  +  \  .  If  the 
convergence  criteria  are  not  met,  cycle  (a)  through  (e)  is  repeated  until  they 
arc  met. 


In  this  report  a  spatially  varying  pseudo-time  step  At  is  used  that  can  be 
interpreted  as  an  attempt  to  use  a  more  uniform  pseudo-Courant  number 
throughout  the  field.'^  A  spacially  varying  At  can  be  effective  for  physical- 
grid  spacings  that  vary  from  very  fine  to  very  coarse  grids,  a  situation  which 
particularly  occurs  for  the  two  Re  =  50  cases.  Pulliam  and  Steger“'  mention 
that  a  spatially  varying  physical-time  step  Al  has  been  used  by  a  number  of 
researchers  to  obtain  solutions  of  steady-state  compressible  lluid  flows. 
Gorski""  applies  a  spatially  varying  pseudo-time  step  for  obtaining  solutions  of 
steady-state  incompressible  fluid  flows.  It  appears  logical,  therefore,  that  in 
this  study,  which  obtains  the  solution  of  a  “steady-state”  incompressible  fluid 
flow  at  each  physical  time  step  At  by  marching  through  pseudo-time  steps  At, 
a  spatially  varying  pseudo-time  step  At  be  used.  This  time  step,  used  in  liq. 
(50),  is 


At  = 


(-56) 
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and  scales  directly  with  the  area  and  aspect  ratio  of  a  physical  grid  cell  ( |7  j  is 
the  area  of  a  cell  in  physical  space).  The  value  of  the  parameter  will  be 
given  later. 

The  form  of  Eq.  (56),  given  by  Hodge, is  borrowed  from  Thompson 
and  Shanks  (Appendix  D)''^  where  it  appears  somewhat  disguised  in  a 
discussion  relating  the  artificial-compressibility  procedure  and  a  Poisson 
equation  for  the  pressure. 

Thompson  and  Shanks,  who  solved  the  time-dependent,  two- 
dimensional  Navier-Stokes  equations  for  the  viscous  fluid  flow  about  a 
hydrofoil  at  a  free  surface,  used  artificial  compressibility  to  obtain  the 
pressure  at  the  free  surface  and  the  hydrofoil,  and  they  employed  a  Poisson 
equation  for  the  pressure  in  the  interior  of  the  flow  region.  In  the  present 
report,  in  contrast,  artificial  compressibility  is  used  to  obtain  the  pressure  in 
the  interior  of  the  flow  region,  and  boundary  conditions  for  the  pressure  are 
used  to  obtain  the  pressure  at  the  free  surface  and  other  boundaries. 

For  the  case  Re  =  50,  Fr  =  0.356  only,  a  free  surface  instability 
developed  at  t  ~  3.52.  Starting  at  this  time,  upon  convergence  at  each  new 
physical  time  (level  n  4-  1),  the  values  for  f  ~  u,  v,  P,  Y  at  the  free  surface 
point  i  are  finally  given  in  terms  of  their  corresponding  converged  values  /: 

L  =  -1-2  +4(/;+,  -i-/;_i)  +  \m  .  (57) 

This  filtering  process  was  developed  by  Shapiro’'*  and  was  used  by 
Longuet-Higgins  and  Cokelet,'^  among  others,  to  eliminate  numerical 
instability  at  the  free  surface.  For  Re  =  50,  Fr  —  0.356,  filtering  was  applied 
at  the  free  surface  from  x  =  1.0  to  2.0  during  the  time  t  =  3.52  to  4.75  and 
from  .V  =  1.0  to  2.94  during  the  time  t  =  4.75  to  9.52. 

The  kinematic  condition  (36)  was  used  for  all  cases  except  Re  =  50, 
Fr  =  1.125  from  /  =  2.40  to  4.06.  During  that  time  Eq.  (37)  was  used  because 
of  the  need  to  crowd  more  free-surface  points  in  the  .^-direction  to  resolve  the 
steep  wave  generated  for  this  case. 

The  relaxation  factors  were  chosen  to  be  1.85  for  Eqs.  (48)  and  (49)  and 
1.0  for  Eqs.  (51)  and  (52).  The  parameters  (cj,  €2)  =  (0.01,  0.01)  were 
applied  to  Eqs.  (48)  and  (49)  and  (rj,  €2)  =  (0.03,  0.01)  to  Eqs.  (51)  and  (52). 

The  adaptive  grid  parameters  {A,  B)  =  (2.0,  2.0)  were  used  for  all 
cases  except  Re  —  50,  Fr  =  1.125  from  t  =  3.10  to  4.06.  During  this  time, 
{A,  B)  were  gradually  increased  to  (3.8,  3.8). 

For  the  two  Re  =  50  cases,  the  parameters  (zV,  ^Tre/)  were  gradually 
increased  from  (10“®,  460.8)  to  (0.0003,  2560)  during  the  time  /  =  0  to  0.019. 
F'or  these  two  cases,  the  parameters  (A/,  Airey)  =  (0.0003,  2560)  did  not 
change  for  the  remainder  of  their  time  histories,  except  that  A/  was  cut  back 
to  0.00015  for  Fr  =  0.356  from  t  —  3.52  to  5.02  and  that  the  parameter  set 
(0.0003,  2560)  for  Fr  =  1.125  was  reduced  to  (0.000025,  1920)  from  t  =  4.0  to 
4.06.  The  total  time  span  extends  to  /  =  9.52  and  t  =  4.06  for  Fr  =  0.356  and 
Fr  =  1.125,  respectively. 


I  or  Re  —  10,  /•>'  -  1.125,  the  spatially  varying  pseudo-time  step  Ar  was 
not  used  because  the  grid  in  physical  space  remained  more  uniform  in  time 
than  it  did  in  the  other  cases,  since  the  How  gradients  are  smaller  for  the  case 
Re  =  10  than  for  the  other  cases,  h'or  Re  =  10  the  parameters  (A/,  Ar)  were 
gradu  ;lly  increased  from  (lO"'.  1.6)  to  (0.00015.  8.0)  during  the  time  /  =  0  to 
0.007.  I'he  parameters  (0.00015,  8.0)  were  then  retained  until  the  end  at 
t  =  9.04. 

I'or  the  two  Re  ~  50  cases  a  maximum  of  three  to  eight  pseudo-time 
step  iterations  per  jihysical  time  step  near  t  ~  i)  were  needed.  Iteratii  ns  b)r 
hr  =  ()..i56  then  gradually  leveled  ol'f  to  i>ne  to  two  iterations  per  physical-time 
step  for  the  remainder  of  the  time.  Iterations  for  hr  =  1.125  leveled  off  to 
two  but  then  increased  to  two  to  live  iterations  per  physical-time  step  as  the 
wave  steepened.  1  he  Re  -  10  ease  rei.|uired  fewer  iterations  i‘)er  physical-time 
step  than  the  Re  ■---  50  cases. 

Mass  conservation  was  monitored  in  the  How  held  for  the  three  How 
eases.  I'he  rletails  are  here  omitted. 

The  computer  time  used  on  the  C'ray-XMP  24  for  Re  —  50,  h'r  =  ()..s56 
and  Re  -  50,  hr  -  1.125.  which  were  run  until  t  —  9.52  ami  4.06. 
respectively,  were  approximately  140  and  65  minutes,  respectively.  Re  =  10, 
which  was  run  until  t  =  9.04.  needed  75  minutes. 


I'o  build  confidence  in  the  numerical  method  used,  the  case  of 
Re  =  50,  hr  =  0..756  was  also  computed  with  a  grid  which  had  four  times  the 
number  t'f  cells  of  the  original  grid,  that  is,  (.713x269)  points.  The  results, 
which  are  discussed  later,  agree  ciuite  well  with  those  of  the  coarser  grid  of 
( 157x135)  |xiints. 

'The  relationship  among  Ar,.,,^,  A/,  the  fineness  of  the  grid,  and  the 
convergence  criteria  still  must  be  explored  for  maximum  efficiency  and 
accuracy  of  the  numerical  method  used  in  this  work. 


An  auxiliary  ciuantity  of  interest  in  the  generation  of  vorticity  at  the  free 
surface  is  the  distribution  of  surface  vorticity  For  a  given  free  surface 
y  =  F(.v.  r)  and  surface  velocity  components  and  v’^,  the  surface  vorticity  is 
given  by  the  well-known  formula  “twice  the  surface  curvature  times  tangential 
velocity”,  or  in  this  report’s  notation 


f  yj)-  ’ 


In  eamiputational  space.  y\,  and  y„  must  be  replaced  by 

y  ^v 
y, 

V',  -  —  -  -  - —7 — 

•Vs  (,Vf)' 


Numerically,  an  alternative  way  of  computing  o.y  can  be  chosen  by  the 
definition  of  .c  -  iiv/iix  -  i)n/i>y  itself.  Since  in  this  case  information  from 
the  interior  points  is  incorporated  into  the  numerical  scheme,  the  result 
appears  more  accurate  Ilian  with  the  use  of  bep  (58)  in  which  y^  aiiiiears,  a 
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quantity  diflicult  to  obtain  with  snt'licient  accuracy.  In  Fig.  5  the  discrepancy 
between  the  two  methods  is  shown  for  the  special  case  Re  =  50,  Fr  =  0.356 
at  t  =  5.02  for  both  the  coarser  and  the  finer  grids.  A  decision  as  to  which 
formula  is  more  accurate  cannot  yet  be  made. 


RKSUf/rS 

Numerical  computations  were  performed  for  Re  =  10  and  50, 
Fr  =  1.125,  and  for  Re  =  50,  Fr  =  0.356  from  the  initial  position 
-Vi,  =  0.5,  Vv  =  ~d  =  —3.  I'he  Fronde  numbers  were  selected  to  coincide  with 
those  of  Sarpkaya  et  al.'^  and  Telste.'^ 

In  the  case  of  Re  =  10,  Fr  =  1.125  diffusion  dominates  convection  so 
much  (the  fluid  is  so  ‘'viscous”)  that  the  maximum  elevation  v  =  0.307  at 
t  =  4.0  caused  by  the  vortex  motion  is  small  compared  to  the  corresponding 
inviscid-llow  case  of  Sarpkaya  et  al.''^  and  'I'elste.^  The  free  surface  returns, 
after  the  maximum  elevation  has  been  reached,  monotonically  to  the  state  at 
rest.  Fig.  6  displays  the  vector  field  of  the  velocity  at  f  =  4.0,  and  F'igs.  7  and 
8  show  equi-  vorticity  lines  at  /  =  4.48  and  t  =  9.04.  The  vorticity  distribution 
on  the  free  surface  is  seen  in  Fig.  9.  The  data  for  the  positions  of  the  vortex 
are  recorded  in  Table  1.  fhe  path  is  plotted  in  Fig.  10  and  compared  with  the 
solution  for  the  potential  flow  with  a  flat  surface,  Eq.  (1).  In  viscous  flow,  the 
center  of  the  vortex  can  be  defined  either  as  the  place  of  extremal  vorticity  (in 
this  case  of  minimum  vorticity)  or  as  the  center  of  the  whirl  (center  of  the 
nested  streamlines).  The  latter  definition,  however,  depends  on  the  choice  of 
the  reference  frame.  In  Fig.  6  the  reference  frame  is  fixed  to  the  undisturbed 
free  surface  and  the  vortex  is  moving  relative  to  this  frame.  If  the  reference 
frame  is  fixed  to  the  vortex  center,  the  position  of  the  center  of  the  whirl 
(moving  relative  to  the  reference  frame)  shifts  closer  to  that  of  the  minimum 
vorticity  (see  next  paragraph  for  an  example).  Even  then,  these  two  locations 
do  not  necessarily  coincide  as  the  analytical  solution  for  a  decaying  vortex 
dipole  demonstrates.”^’  Fig.  10  shows  that  the  path  of  the  point  of  minimum 
vorticity  is  closer  to  the  axis  of  symmetry  and  closer  to  the  free  surface  than 
the  path  of  the  whirl’s  center.  The  phenomenon  of  “rebounding”  is  observed 
in  both  cases,  that  is,  the  turning  away  of  the  vortex  from  the  free  surface. 
'Fhe  elevated  surface  returns  to  the  state  at  rest  without  oscillation.  At  the  last 
computed  time  t  =  9.04,  |  |  is  0.2383,  diminished  from  the  initial  value 

I  ^’min  1  =  39.43.  In  Fig.  1 1  |  I  is  plotted  against  time.  The  curve  follows 
closely  the  l/t-decay  according  to  the  Hamel-Oseen  solution,  Eq.  (25), 
indicating  that  numerical  diffusion  is  minimal. 

I'or  Re  =  50,  Fr  —  1.125  the  situation  is  quite  different.  The  surface 
elevations  at  three  different  times  in  F'igs.  12  through  16  reveal  a  much 
stronger  effect  of  the  ascending  vortex,  and  a  local  depression  of  the  surface, 
called  a  “scar”  by  Sarpkaya  and  Henderson,®  is  now  visible  that  was  not 
apparent  in  the  case  of  Re  =  10,  Fr  =  1.125.  The  surface  elevations  for 
Re  =  50,  Fr  =  1.125  are  closer  to  the  curves  for  the  inviscid  fluid.  In  fact. 
Figs.  12,  15,  and  16  may  be  compared  with  F'igs.  5a,  5c,  and  5e  of  Sarpkaya  et 
al.'^  The  dimensionless  time  T  in  their  work  is  related  to  this  report’s  time  t  by 
T  ~  t  —  3.  Although  it  takes  a  little  longer  for  the  viscous  fluid  flow  to  reach 
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states  comparable  to  those  of  the  iiiviscid  Iluid,  the  qualitative  agreement 
between  the  two  cases  is  good  with  regard  to  free-surface  height  at  the 
centerline  and  the  locations  of  whirl  center  and  scar.  The  same  statement  can 
be  made  about  'I’elste’s  data.'"*  'I’hey  can  be  compared  with  those  in  this 
report,  if  one  considers  the  relation  between  the  dimensionless  time  and 

t  to  be  —  2j(f  +  2).  I'he  I'roude  numbers  are  related  by  Fr-fehte  =  ^nFr. 

In  I'igs.  13  and  14  a  section  of  the  How  field  of  Fig.  12  is  compared  for 
two  reference  frames.  I’he  reference  frame  fixed  to  the  center  of  the  whirl 
(Fig.  14)  shows  that  the  center  is  shifted  to  the  left,  but  the  new  center  docs 
not  reach  the  location  of  \  v.',,,;,,  j .  'fhe  eciui-vorticity  lines  are  displayed  in  l-'igs. 
17  through  19.  'fhe  merger  of  the  vorticity  field  of  the  vortex  with  that  at  the 
free  surface  is  clearly  visible.  The  vorticity  distribution  on  the  free  surface  is 
given  in  Fig.  20.  The  path  of  the  vortex  center  is  plotted  in  Fig.  21  and  the 
positions  of  the  vortex  center  put  together  in  Table  2.  The  movement  of  the 
vortex  is  straight  up  and  no  deflection  due  to  tiie  free  surface  is  observed.  'I’he 
decrease  of  |  1  with  time  is  given  in  Fig.  22  and  follows  the  l//-law.  The 

computations  were  stopped  at  /  =  4,06  when  convergence  thereafter  could  not 
be  obtained  for  the  grid  used.  However,  a  negative  horizontal  velocity 
component  at  the  free  surface  indicates  that  the  trend  toward  a  constriction  is 
present,  as  it  is  in  the  case  of  the  inviscid  Iluid.  The  mounded  shape  of  the 
free  surface  around  the  symmetry  line  was  still  rising  when  the  computation 
was  terminated. 

A  comparison  of  Fig.  20  with  I'ig.  22  shows  that  the  surface  vorticity 
becomes  larger  than  | -Cn,!,,  |  of  the  vortex  from  approximately  I  =  4  on.  This 
is  an  interesting  situation  because  it  means  that  the  vorticity  concentration  of 
the  tlow-gcnerating  vortices  is  less  than  the  surface  vorticity  generated  by  the 
vortices.  'I'he  statement  contrasts  the  result  obtained  for  the  flat  slip-surface 
by  Peace  and  Riley, that  the  flat  surface  is  always  a  sink  of  positive  vorticity 
(or  a  source  of  negative  vorticity)  because  =  0. 

The  final  case  computed  was  Re  —  50,  Fr  =  0.356.  The  lower  Fronde 
number  means  less  disturbanee  of  the  free  surface  by  the  vortex  motion.  This 
smaller  disturbance  is  observed  in  Figs.  23  through  25.  On  the  centerline  the 
free  surface  reaches  a  maximum  elevation,  then  falls  to  a  minimum,  climbs 
again  to  a  maximum,  and  comes  to  the  state  at  rest.  In  other  words,  the  free 
surface  oscillates.  This  oscillation  is  in  contrast  to  the  case  of 
Re  =  10,  Fr  =  1.125  that  showed  only  an  up  and  down  movement  of  the  free 
surface  at  the  center  line.  The  scar  is  now  pronounced,  with  the  high  surface 
curvature  that  results  in  high  surface  vorticity  (Figs.  26  through  31).  From 
this  concentration  of  vorticity  a  secondary  vortex  develops  which  becomes 
visible  at  /  =  5.02  in  front  of  the  primary  vortex  (Fig.  28).  This  figure  was 
also  computed  with  the  liner  (313x269)  grid  in  I"ig.  29.  The  two  figures  agree 
well  except  for  the  low-level  positive  values.  A  little  later,  at  t  =  6.52,  the 
secondary  vortex  has  placed  itself  directly  in  front  of  the  primary  vortex 
farther  away  from  the  free  surface  (I'ig.  30).  The  computations  were 
continued  without  any  numerical  difficulties  and  were  stopped  at  /  =  9.52.  In 
I'ig.  32  the  free-surface  vorticity  is  displayed,  and  in  I'ig.  33  the  path  of  the 
vortex  center.  In  'I’able  3  the  position  data  are  recorded.  In  contrast  to  the 
previous  case  of  Re  —  50,  Fr  -  1.125,  the  path  of  the  vortex  center  now 
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shows  a  rebound  from  the  free  surface.  The  decrease  i>f  |  I  with  time  is 
siiowii  in  J'ig.  .^4.  I  hc  agreement  between  the  curves  in  lugs.  22  and  34  is  very 
good,  and  the  1 //-law  is  observed.  I'his  means  that  the  decay  of  tlie  vortex- 
center  vorticity  is  not  intluenced,  or  barely  so,  by  the  presence  of  the  free 
surface.  Here  again,  the  ilata  are  compared  willi  those  for  the  liner  grid,  and 
almost  no  distinction  can  be  observed. 

I'he  velocity  field  in  I'ig.  2,"^  can  be  compared  with  that  for  an  inviseid 
Iluid  as  shown  by  I'elste.'^  'I’he  scars  are  less  pointed  in  the  case  of  a  viscinis 
lluid  than  in  the  case  of  an  inviseid  one,  and  the  position  of  the  vortex  center 
is  knver  in  the  viscous-tlow  case  than  in  the  inviscid-lluid  one. 

The  intlueirce  of  surface  tension,  the  case  of  an  oblique  approach  of  a 
vortex  pair  toward  the  free  surface,  and  an  attempt  to  compute  Hows  with 
higher  Reynolds  number  will  be  treated  in  a  forthcoming  paper. 
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Table  1 


Pesilions  .v.  of  i i,, !  aiul  center  of  whirl  as  a  Innction  of  lime  /  for 

AV  -  11).  1.125. 


w  liirl 


^  j ,  0 

0 . 50000 

-  3 . 00000 

0 . 53Z40 

-2.87100 

0 . 52 

0 , 53'?50 

-2 . 55443 

0.72 

-  2 . 6 1 

1 . 00 

0 . 59450 

-2 . 23920 

0 . 04 

-2.37 

i  .  4c 

0.64670 

-2 . 07500 

0 . 95 

""  Z  •  iL 

2 .  OZ 

0 . 74830 

-1.91 060 

1  .08 

-2 . 02 

Z .  5' 

0 . 79907 

-i.7997£ 

1  .17 

- 1 . 9- 

04 

0 . 84954 

- 1  .  /-.8650 

t  .24 

- 1  .  :36 

•.  S 

0 . 9.5000 

-1.63019 

1 

•1.82 

4 . 00 

1 .00017 

-1  .52198 

1  . 2:9 

-1  .■'9 

4 , 43 

1.05031 

- 1 . 52026 

1 .44 

- 1  . 84 

5 .  '"'Z 

1 . 15026 

-1.47129 

1  , 2'0 

-  1  .  '2.-. 

*=  ..  50' 

1 . 2002 1 

-1 .4/590 

1 . 54 

- 1  . 

■  .  04 

1 .25013 

- 1 . 43300 

1 . 5:3 

-  1 . 9  3 

1 . 34996 

-1.44439 

; .  59 

-1.95 

7 .  C»(!* 

1 . 39995 

- 1  .50670 

1  . 63 

7.54 

1 . 44993 

- 1 . 52009 

1 . 70 

-2.‘j4 

3 . 0  2 

1.49997 

-1.58206 

1.75 

-2.14 

8 . 50 

1 . 54994 

-1.59290 

1 , 78 

-2.  \2 

V .  04 

1 . 60000 

- 1 .65366 

1.79 

-2.16 
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Table  2 


Positions  .v,  v  of  |  [  and  center  of  whirl  as  a  function  of  time  t  for 

Re  =  50.  /•>-  ’=  1.125. 


U  lill’ 


t 

X 

y 

X 

y 

0 ,  (j 

0 . 5000 

-3.0000 

0 . 007 

0 . 5079 

-2-9920 

0 . 0 1  '5 

0 . 505'S 

-2.9861 

0 .  '04<j 

0 . 6043 

-2.9616 

0. 160 

0.5077 

-2.8419 

0.5i"i) 

0.5159 

-2.4853 

0 . 55 

-2.50 

1 . 000 

0 . 5Z35 

-2.0527 

0 . 58 

-2.06 

1  . 4;-;0 

0 . 5Z62 

- 1 . 6- 1 96 

0,62 

- 1 . 64 

2.01:0 

0 . 5ZS0 

-1.1672 

0 . 65 

-1.19 

2.401 

0 . 5604 

-0.8926 

0 . 66 

-0.91 

2 . 500 

0 . 5'367 

-0.3115 

0 . 68 

-0.S3 

.  040 

0 . 5558 

-0.4561 

0.70 

-0.49 

0 . 1 00 

0.5705 

-0.4207 

0.71 

-0.46 

3 . 520 

0 . 5677 

-0.2185 

0.70 

-0.23 

3.820 

0.5751 

-0 . 0803 

0 . 68 

-0.13 

:]! .  940 

-0.08 

4 . 000 

0.5798 

-0.0249 

0 . 68 

-0.05 

4 . 050 

0.5890 

-0.0284 

4 . 060 

0 . 5 '^09 

-0.0265 

u  *  <66 

-0.04 
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I'ahic  3 


Posiiioiis  A,  V  of  (v-niiiij  and  center  of  whirl  as  a  I'uiictiou  til'  time  /  for 
Ri>  -  51),  /•>  -  0.356. 


V 

y 

y 

0 . 0 

0 .  ^I'OOO 

“3 . 0000 

0 . 00  7 

0 . 5079 

-2.9920 

0.019 

0 . 5056 

-2.9861 

0 . 040 

0 . 504S 

-2 „ 96 16 

0 . 1 60 

0.5077 

-2.8419 

' ) . 5Z0 

0.5159 

-2.4854 

0 . 55 

-2.50 

1 . 000 

0.5238 

-  2 . 05'5  7 

0 .  f<8 

-2 . 06 

1 . 4:00 

0 . 5272 

-1.6353 

0 . 6Z 

-1.65 

2.000 

0 . 5608 

-1.2605 

0 . 66 

- 1  . 28 

2 . 500 

0 . 6255 

-0.9778 

0.71 

-1  .02: 

s .  040 

0.7208 

-0 . 7876 

0 . 80 

-0.88 

;! ,  5?0 

0 . 9058 

-  (j  .  7 1 60 

0.92 

-0.:-::5 

0 . 760 

0.9999 

-0.  7:288 

1 .01 

-0.:?5 

4 . 000 

1 . 1247 

-0.7406 

1 . 1 0 

-0.88 

4.240 

1 .2189 

-0 . 7  796 

1  . 20 

-0.'''0 

4 . 400 

1 . 34:39 

-0.7836 

1.32 

-0.92 

4.750 

1 .4:579 

-0.8160 

1 . 42 

-0.93 

5 . 020 

1.5314 

-0 . 8:454 

1  .50 

-0.94 

5 . 500 

1 . 6566 

-0 . 8678 

1 . 64 

-0 . 96 

rj  .  040 

1.7811 

-0 . 8855 

1 . 75 

-0.46 

6 . 5,20 

1.8752 

-0.8716 

1 . 84 

-0.98 

7 .  ( 'Oi'i 

1 . 9690 

-0 . 8950 

1.94 

- 1 . 00 

7 . 540 

2 . 0942 

-0 . 8:90:5 

c-  ■ 

- 1 . 00 

0 .020 

2.1881 

-0.8924 

2.14 

- 1 .  Oi, 

.  500 

Z . ZS20 

-0.8995 

2.24 

-1  .''2 

4.040 

2.3760 

-0.9107 

Z  a  i  :'' 

- 1 . 0-1 

0 , 520 

2.4701 

-0.9211 

2 . 43 

- 1  . 63 
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Fig.  1.  (  cl)  Sketch  of  the  How  situation,  (b)  Mapping  of  the  physical  s 
onto  the  computational  space,  (c)  “I-’our-color”  scheme. 
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Vector  lickl  of  the  velocity  for  Re  ~  10,  Fr  =  1.125  ;it  / 


0  1 


2  3  4  5 

X 


Fig.  8.  Equi-vorticity  lines  for  Re  =  10,  Fr  =  1.125  at  t  =  9.04.  I'hc 
of  the  contours  are  —0.15,  —0.1.3,  etc.  Solid  lines  represent  ne 
clashed  lines  positive  data. 
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Fig,  11.  Decrease  dI  [ ]  willi  time  I'or  Re  =  10,  I'r  ~  1.125. 


inch  -  5  ussr  units 


- - - T'TV'/ 

/ // ^ 
j' 

///  / 

y/*y*/‘  /  /^  /  ^  ^  ^  ^  ^  ^  y,  y 

'’y'y/ / /-  y  ^ 

^  y'  y  ^  y  ^  ^  y'  ^  y  ^ y  y  y  y  yy 
^  '^y^'^y^yyyyyyyyyy 

y'  y  y y  y  y  y  y  y  y  y . 
''''''''/-'yy/'yyyyyyyy^ 
y  ^  y^  yyyyyyyy^ 
^  yy  y  yyyyyyyy 

////  //////  yyy  yyyyyy^^ 

^y'^y  yyyyyy^^ 
////  f  {''''''  y  ^  y /y  yyyyy^yy, 

/ff f  /  -'  '^'^'^/yyy yyy yyyyy, 

////  /  /  /  /  /  /  /  /  / /yyyyyyy^ 
III!  I  /  /  //  //  //  //  /  yyy  yyy 

'v7///////////r;^;^/^;;y 


///  /  /  /  / 
J  !  !  I  t  i  / 
>1111/1. 
'I  J  I  I  J  J 


I 


\''  //i'  i  III i n  ( ( i '' ''' 

>  ^  ^  0  0  7  \''  \\  V  \\  \\  V  V  V  i  0  /  /  '  '  '  '  •■ 

'  \  “v  “v  \^\V\\^\\\V^\VNyKV7>A^ — ' 

M 


->''S\\\\\\ 

V  ^  \  \  ,\  ,\ 

" '  II 1 1 1 1 


y^l^ilhhlh 


13.  Scclioii  ol  lli^  ll^’w 


Same  How  Held  as  in  I'i”.  13  but  with  rclerenee  I'rame  fixed  to  the  vortex  eenter. 


Fig.  16.  Vccloi  liclii  ol'  the  velocity  lor  /vV  50,  />'  ^  i,125  at  t  4.06. 


Fig.  20.  I'lcc-surl'acc  vorticily  l\>r  Re  50,  Fr  =  1.125  at  (a)  t  —  2.50,  (b)  t  =  3.52,  and  (c)  t  =  4.00. 


of  the  vortex  eenter  for  Re  =  50,  Fr  =  1.125,  defined  (a)  by  ]  j ,  (b)  by  the  center  of  the  whirl, 
oinparison,  l.amb’s  .solution,  Fq.  (1),  is  added.  Also,  the  free  surface  at  the  last  computed  time 
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Fig,  22.  Dccrcjisc  of  i !  witli  lime  for  Re  =  50,  /•'/■  =  1.125. 
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Fig.  23.  Vector  (icltl  of  llic  velocity  lor  Re  =  50,  Fr  —  0.356  at  /  =-  3.52 
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51 


Fig.  30.  Equi-vorticity  lines  for  Re  =  50,  Fr  =  0.356  at  t  =  6.52.  The  values  of  the  contours  are 
±0.3,  ±0.9,  ±1.5,  etc.  Solid  lines  represent  negative,  dashed  lines  positive  data. 


53 


Fig.  32.  Frcc-siirt'ace  vorticity  for  Re  =  50,  Fr  =  0.356  at  (a)  t  =  3.52,  (b)  t  =  5.02,  and  (c)  /  =  6.52. 


;  vortex  center  lor  Re  =  50,  I  r  ~  0.356,  defined  (a)  by  |  |  for  bolb  fme  and  coarse  grids, 

center  of  the  whirl,  for  comparison.  Lamb’s  solution,  Eq.  (1),  is  added.  The  two  curves  (a) 
indistinguishable. 
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